Molecular genetics and quantitative traits divergence among populations of Eothenomys miletus from Hengduan Mountain region

Abstract An important objective of evolutionary biology has always been to grasp the evolutionary and genetic processes that contribute to speciation. The present work provides the first detailed account of the genetic and physiological adaptation to changing environmental temperatures as well as the reasons causing intraspecific divergence in the Eothenomys miletus from the Hengduan Mountain (HM) region, one of the biodiversity hotspots. One hundred sixty‐one E. miletus individuals from five populations in the HM region had their reduced‐representation genome sequenced, and one additional individual from each community had their genomes resequenced. We then characterized the genetic diversity and population structure of each population and compared the phenotypic divergence in traits using neutral molecular markers. We detected significant phenotypic and genetic alterations in E. miletus from the HM region that were related to naturally occurring diverse habitats by combining morphometrics and genomic techniques. There was asymmetric gene flow among the E. miletus populations, indicating that five E. miletus populations exhibit an isolation‐by‐island model, and this was supported by the correlation between F ST and geographic distance. Finally, P ST estimated by phenotypic measures of most wild traits were higher than differentiation at neutral molecular markers, indicating directional natural selection favoring different phenotypes in different populations must have been involved to achieve this much differentiation. Our findings give information on the demographic history of E. miletus, new insights into their evolution and adaptability, and literature for studies of a similar nature on other wild small mammals from the HM region.


| INTRODUC TI ON
Early flora and fauna cradleland and refuges are hotspots for biodiversity. Some biodiversity hotspots serve as "evolutionary forewords" that spur fast divergence of tropical plant groupings and the junction of long-distance species distribution, having a significant impact on the establishment and evolution of the world's flora and fauna.
Due to Hengduan Mountain (HM) region has the characteristics with high northwest and low southeast, and its climate, which is characterized by a modest yearly temperature difference and a huge daily temperature difference, a wider range of animal species can survive there (Ren, Jia, et al., 2020). Animals display a variety of phenotypic alterations as a result of selection forces acting on heritable features as a result of geographical and temporal heterogeneity (Leinonen et al., 2008). Animals may go through these phenotypic changes to better fit their environment at the physiological, behavioral, and especially morphological levels. Although phenotypic plasticity has been extensively studied and its significance in adaptation and evolution has been well-discussed, the basic driving mechanisms are still unknown (Kelly et al., 2012;Sommer, 2020).
Comparative analyses of quantitative genetic and neutral marker differentiation have given researchers a way to assess the relative contributions of stochastic genetic drift and natural selection to the explanation of among-population divergence (Leinonen et al., 2008).
In several species, the comparison of quantitative trait across populations (Q ST ) and differentiation at neutral molecular markers (F ST ), commonly referred to as the Q ST -F ST comparison, revealed that natural selection played a significant role in the cause of differentiation in quantitative traits. In several cases, putative F ST and Q ST differentiation in various populations is compared in order to evaluate their evolutionary signatures and discover potential features implicated in local adaptation.
However, raising animals from various populations in a common environment is typically required for estimating the additive genetic variances needed for Q ST (Brommer, 2011;Leinonen et al., 2008). As a result, for some wild species, particularly endangered species, the breeding test for estimating the Q ST becomes impractical. Currently, most species substitute quantitative trait analysis (Q ST ) with phenotypic divergence in a trait (P ST ), and P ST counts are based on phenotypic assessments of a wild trait of several individuals across numerous populations (Brommer, 2011). P ST -F ST comparisons, on the other hand, rely on the unrealistic presumption that nonadditive genetic effects and environmental effects may be reduced and that phenotypic variation equals additive genetic variance (Wójcik et al., 2006).
Eothenomys of subfamily Arvicolinae, which belong to the family Cricetidae in Rodentia, is wildly distributed throughout the Holarctic realm and parts of the Oriental realm (Luo et al., 2004). Longstanding controversy surrounds the precise phylogenetic position of Eothenomys. There are obvious differences in morphological and molecular evolution of the Eothenomys (Liu, 2007;Liu et al., 2012).
Recently, according to research on the species of Eothenomys utilizing molecular and morphological evidence revealed that Eothenomys has three subgenera, which includes Eothenomys, Anteliomys, and Ermites (Liu et al., 2012;Zeng et al., 2013). It is possible that the dispersal of Eothenomys started late and tended from south to north and west to east during the whole evolutionary process (Liu, 2007;Ye et al., 2002).
E. miletus is a naturally occurring species in HM region (Ren, Liu, et al., 2020;Zhu et al., 2010), and is listed in International Union for Conservation of Nature. E. miletus is one of the representative species for studying the evolution of biodiversity in HM region (Zhu et al., 2008). Our previous studies have shown that there were significant difference in phenotypic morphology (Zhu et al., 2008(Zhu et al., , 2010(Zhu et al., , 2011(Zhu et al., , 2014, energy metabolism, level of serum leptin and expression of hypothalamic neuropeptides (Ren, Liu, et al., 2020), and microsatellites (Zhu et al., 2013) of five E. miletus populations, but our understanding of evolution and adaptation within E. miletus populations is limited due to the lack of genomic studies. In the present study, we use simple genome sequencing and resequencing techniques to explore the degree of genetic differentiation and genetic structure among five E. miletus populations from HM region, as well as compare the quantitation of the P ST based on the collected the morphological data with F ST estimated using sequencing to assess the relative effect of different evolutionary mechanisms on the phenotypic differences among E. miletus populations in HM region.
Finally, we provide literature for the similar studies on other wild small mammals from HM region.

| Subjects and experimental design
From November 2018 to January 2019, the voles (E. miletus) used in this study were caught in five sites with gradually varying altitudes: Deqin (DQ, n = 33); Xianggelila (XGLL, n = 33); Lijiang (LJ, n = 34); Jianchuan (JC, n = 33); and Ailaoshan (ALS, n = 33). Detailed sampling information is provided in the Figure 1a and Table 1. Latitude, elevation, and annual mean temperatures in the study were obtained from local weather bureaus. Animals caught in the wild were weighted immediately killed by anesthesia, and their liver were immediately dissected and frozen in liquid nitrogen. Samples were stored in dry ice and transported to the laboratory of Yunnan Normal University and stored at −80°C refrigerator until assayed. Using the phenol/ chloroform method, the total genomic DNA of the animals was extracted from tissue samples. For each individual, 1-3 μg of DNA was sheared into 200-500 bp fragments with the Covaris system (Gene

| Sample sequencing, read mapping and quality control
We selected the reference genome of Microtus ochrogasternas the prediction of electron enzyme digestion. DNA samples were digested with Rsal (Beijing Baimai Biotechnology Co. Ltd), and the sequence of fragment lengths 464-494 was defined as a SLAFtag. The principles of the selective digestion plan are as follows: (1) The proportion of digestion SLAF tags located in repeat sequences should be as low as possible; (2) the SLAF tags were evenly distributed on the genome as far as possible; (3)

| SNP calling and filtering
In order to estimate the sequencing quality value Q, the reads considered to be of low quality were those with joint and 50% bases with a Q10 value. Q = −10*log10e. Polymorphic SLAF tags were employed for single nucleotide polymorphism locus (SNP) calling in the GATK 3.3.0 (McKenna et al., 2010) and SAMtools v1.2  routines. The SNPs were further filtered to exclude SNPs present in <80% of individuals and those with a minimum allele frequency (MAF) < 5% using PLINK v1.90 (Purcell et al., 2007).
Clean paired-end reads from resequenced individuals were aligned to the resequenced assembled vole reference genome using BWA 0.7.12-r1039. Then, SNPs were identified using GATK 3.3.0, and the clean SNPs were aligned using GATA 3.3.0 hard filter with the following filter parameters: QD 2.0, FS >60.0, MQ 40.0, ReadPosRankSum −8.0, and MQRankSum −12.5. Only SNPs with high credibility were retained for further analysis after the SNPs were filtered by MAF = 0.06 and maximum missing rate = 25%.

| Population structure
Population structure analysis was performed using the ADMIXTURE  (Zhang et al., 2022).

| Genetic diversity
The expected heterozygosity (H e ) and observed heterozygosity (H o ) were calculated using PLINK 1.9 (Purcell et al., 2007) to test the genetic diversity indices of five populations based on highconsistency SNPs, and the Nei diversity index, and polymorphism information content (PIC) were calculated using a customized Perl script (Beijing Baimai Biotechnology Co. Ltd). We used SPSS 26.0 to calculate Pearson's correlation coefficient (r 2 ) between each pair of variables in order to further quantify the impact of environmental variables, such as altitude, temperature, and latitude, on genetic diversity at five geographic populations (Qu et al., 2014).
All environmental variable using in the analysis were downloaded from National Meteorological Science Data Center (http://data. cma.cn/).

| Demographic history reconstruction and gene flow
The maximum likelihood method and a Bayesian statistical model were

| Neutral genetic differentiation and phenotypic differentiation
SNPRelate package in R 3.6.3 was used to calculate pairwise F ST (Zheng et al., 2012), and Prism 9 was used to build a heat map of the pairwise F ST value. Based on Pearson's product-moment correlation, the Mantel test of matrix correspondence (Mantel, 1967) was applied to test correlations between geographic distance, environment distance, altitude distance, temperature distance, precipitation distance, pairwise F ST , and F ST /(1 − F ST ) in order to test the effects of geographic distance and environmental differences on genetic differentiation. This was done using the Ecodist package in R 3.6.3 (Rousset, 1997). On topographic maps of the study area, point-to-point geographic distances were calculated (Browne & Ferree, 2007). Moreover, we gathered environmental data from WorldClim 2.0 for sampling locations using 19 common bioclimatic variables (Fick & Hijmans, 2017). ArcMap 10.2 was used to convert the data (Berhanu et al., 2022). The 19 standard bioclimatic variables that correlate with temperature were utilized as temperature data, while the 19 typical bioclimatic variables that relate to precipitation were used as precipitation data.
To calculate the distance in environment, temperature, and precipitation, we employed the Pearson distance measurement method. General linear regression analysis in R 3.6.3 was used to investigate the relationship between geographic distance and envi- as the response, the geographic distance as the predictor, and the environmental distance as the condition factor to assess isolation by distance (IBD). The geographic distance was utilized as the condition element to investigate isolation by environment, isolation by temperature (IBT), and isolation by precipitation (IBP). Moreover, the distance in altitude between paired sampling sites was calculated.
We created small mammal skull specimens using the Tenebrio molitor larval method (Ren, Jia, Zhang, Wang, & Zhu, 2023). The analysis of the fractured skull specimens was not carried out. One hundred twelve complete skull specimens were maintained at animal speci- Using the SPSS 26.0 program, the body mass and 21 exterior and cranial character data were evaluated. One-way analysis of variance (ANOVA) and LSD post-hoc tests were used to assess group differences in attributes; p < .05 was deemed statistically significant, while p < .01 was deemed extremely significant. Prism 9 was used to perform Highcharts and Boxplot. Using the online Heatmapper, a cluster analysis plot and correlation matrix map between physical characteristics and environmental factors were created.
Divergence at phenotypic traits will be greater than that seen for neutral loci under divergent selection.  (Spitze, 1993), which quantifies the proportion of among-population genetic variance in quantitative traits. We calculate the P ST based on the following equation using EXCEL. (Raeymaekers et al., 2007).
where σ 2 B is the variance between populations, σ 2 W is the variance within populations, and h 2 the heritability. The scalar c expresses the proportion of the total variance that is presumed to be because of additive genetic effects across populations, assuming that environmental variance among samples is randomly distributed or absent and that heritability (h 2 ) within samples is 0.5. The consequences of departure from these assumptions are considered below in the Discussion sections. Phenotypic variance components were estimated following Sokal & Rohlf, 1995

| Sample sequencing, read mapping, quality control, SNP calling, and filtering
We sampled 161 individuals of five E. miletus populations from the HM regions ( Figure 1a; Table 1). After quality control, 363.16 million pair-end reads with an average of 92.23% Q30 and 42.09% GC were produced after quality control (Table S1). The 847,420 SLAF tags of 161 individuals with 9.19× depth coverage, which included 470,440 polymorphism labels, were collected (Table S2). SNP information of each sample is shown in Table S3.
In the whole-genome resequencing experiment of five individuals, we obtained 0.65 Gb of data that translated into mean genome coverage of each individual between 38× and 42×. Moreover, we obtained clean data with average Q20 and Q30 values of 97.72% and 92.85%, respectively (Table S4), and 108,005,364 SNPs were gathered by comparing with the first 40 scaffolds of the reference genome (Table S5).  (Figure 1c; Figure S1). The five groups spread over these locations showed varying degrees of mixed ancestry as K climbed from 5 to 10.

| Population structure
We next clustered individuals using phylogenetic reconstruction ( Figure 1d). The results of this analysis revealed that the grouping of populations, which showed four clusters, JC and ALS population formed a cluster respectively, and part XGLL individuals with DQ individuals formed one cluster, the remainder XGLL individuals with LJ individuals formed one cluster. This is consistent with the observation from our structure analysis and PCA.

| Demographic history and gene flow
To estimate the pairwise relative migration rates and direction between pairwise populations, we employed the Migrate-N (Figure 2a). and 30 Ka years ago, the two period of low temperature in history (Howard, 1997).  Blunier, 2001). After the periods of fluctuation, the N e decreased. We also found that ALS is diverging from the others no later than 10 5 years ago while the others are together until very recently, which is consistent with high geneflow between DQ, XGLL, and LJ. But a bit surprising for JC looking divergent on PCA, NJ tree, and F ST , which cause by diverging recently.

| Neural genetic differentiation and phenotypic differentiation
The There were extremely significant differences in body mass as well as twenty external and cranial characters, expect AVL, between five populations (p < .01; Figure 4A  between most phenotypic traits and environment factors, which had positive correlation with annual environment temperature, and had negative relationship with altitude and latitude (p < .05; Figure 4D).
We further calculated the pairwise P ST of all phenotypic traits between five populations, and compared with the pairwise F ST .
First the results of violin diagram show that the probability of P ST more than F ST is large (Figure 5a). Moreover, the results of independent sample t-test showed that P ST of all tested traits was significantly greater than F ST (p < .01) . From the two-way clustering heat map of P ST /F ST value, several interesting findings have emerged. First, most of pairwise P ST of phenotypic traits were higher than the pairwise F ST (Figure 5b;  Table 3).

F I G U R E 4
Group differences in body mass (a) and 21 phenotypic traits (b) of five E. miletus populations from HM region. Data were analyzed by one-way ANOVA followed by the LSD post-hoc test. Significant differences were indicated by different alphabetic letters. One-way clustering heat map based on the body and skull traits in E. miletus (c). The correlation matrix between altitude, annual average temperature, and latitude with 22 phenotypic traits (d). ALS, AiLaoShan; BM, Body mass; LTRL, lower tooth row length; AVL, auditory vesicle length; BL, body length; CBL, cranial basal length; CCL, covering cap length; CD, chest depth; CH, cranial height; CL, cranial length; CW, chest width; DQ, DeQin; EL, ear length; ESL, eye socket length; EW, ear width; FLL, fore limb length; HLL, hind limb length; IB, interorbital breadth; JC, JianChuan; LJ, LiJing; NW, neurocranium width; PNL, pillow nose length; T 1 L, tail length; T 2 L, torso length; UTRL, upper tooth row length; XGLL, XiangGeLiLa; ZB, zygomatic breadth.

| DISCUSS ION
Phenotypic changes at the morphological, physiological, and behavioral levels to adapt the diverse environment in HM region were found in E. miletus (Ren, Liu, et al., 2020;Zhang et al., 2019;Zhu et al., 2014). Genetic variations were also found in five E. miletus populations from HM region in this study, and although sharing a similar demographic history, the populations had a clear genetic structure.
According to the result of population structure, there were four clusters in genetic level, which grouped together a part of Xianggelila individuals and Jianchuan population, and the remainder of Xianggelila individuals and Deqin population, and Jianchuan population as well as Ailaoshan population formed a single cluster, respectively. These are different from the statistic of phenotypic variations, which clustered together the Deqin population and Xianggelila population, and grouped together the Lijiang population, Jianchuan population, and Ailaoshan population (Ren, Jia, et al., 2020;Zhang et al., 2019).
High genetic variation can serve as the basis for adaptability to environmental change through natural selection, which is essential to the long-term survival of populations (Bijak et al., 2018;Ellegren & Galtier, 2016), as seen in this study with E. miletus. Geographical   (Reese et al., 2001) and southern Appalachia (Browne & Ferree, 2007). The isolation-by-island model predicts that there was no relationship between the distance separating populations and the amount gene flow, in contrast to the isolation-bydistance mode, which asserts that populations separated by shorter distances will experience higher rates of gene flow than populations separated by longer distances (Browne & Ferree, 2007). Isolationby-island concept typically manifests in animals whose habitat is cut off by an extreme environment, and in those species, the distributions of the sub-populations are typically entirely discontinuous in that environment (Qu et al., 2004). distance in the present study (Browne & Ferree, 2007). Phenotypic changes at the morphological levels to adapt the diverse environment in HM region were also found in E. miletus in this study. This is consist with the previous studies (Ren, Jia, et al., 2020;Zhang et al., 2019).
The clustering of populations by genetic data differed from that by phenotypic data, which may be the result of long-term adaptation of different phenotypes to different habitats (Liu, 2007;Liu et al., 2012).
Moreover, morphological changes had negative correlation with altitude and latitude, and positive correlation with annual environment temperature, indicating that morphological traits of E. miletus dose not obey the Bergmann's rule (Ashton et al., 2000;Bergmann, 1847).
No data were available to estimate the genetic variances of traits in this study due to the fact that the animals in this study are wild, but we can determine the effect on P ST under different h 2 conditions. We further calculated the P ST value using four different heritability estimates (0.25, 0.5, 0.75, and 1), based on the assumption that there is no environmental variance. The graphs in Figure 6 showed the value that  Table S7. These conditions are unlikely to be compatible in nature because nonheritable variance should be environmentally pliable (Wójcik et al., 2006).

| CON CLUS ION
In this study, we investigated the widely dispersed E. miletus in the HM region and used population genomic techniques to provide insights on its differentiation, adaptation, and history. In conclusion, our data show that E. miletus from the HM region exhibits phenotypic and genetic alterations related to naturally occurring diverse environments. It is interesting to note that there are differences between phenotypic clusters and genetic change patterns. Furthermore, phenotypic and genetic changes are linked to environmental factors, such as latitude, altitude, and average annual temperature, and phenotypic traits are more influenced by environmental factors; however, it is still unknown whether other environmental factors may also have an impact on phenotypic and genetic changes. Additionally, the significant biological stratification brought on by the tectonic uplift of the HM region during the late Pliocene results in spectacular topography, which has an impact on the asymmetric gene flow patterns found in E. miletus. Five E. miletus populations demonstrate an isolation-by-island model, which is supported by gene flow and a link between F ST and geographic distance. Last but not least, P ST estimates for the majority of wild traits are higher than differentiation at neutral molecular markers, indicating that directional natural selection favoring various phenotypes in various populations was likely involved in achieving this much divergence. Our findings provide as a foundation for studies on other HM region wild small animals.

ACK N OWLED G M ENTS
We appreciate the assistance of all the Physiological Ecology Group members at Yunnan Normal University in carrying out the experiments and discussing the findings.

FU N D I N G I N FO R M ATI O N
The Yunnan Ten Thousand Talents  Scientific Research (SXBYKY2022123) support for this work.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no competing financial interests.

O PE N R E S E A RCH BA D G E S
This article has earned an Open Data badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available at https://doi.org/10.5061/ dryad.kkwh7 0s8b.